Impact of basic network motifs on the collective response to perturbations

Many collective phenomena such as epidemic spreading and cascading failures in socioeconomic systems on networks are caused by perturbations of the dynamics. How perturbations propagate through networks, impact and disrupt their functions may depend on the network, the type and location of the perturbation as well as the spreading dynamics. Previous work has analyzed the retardation effects of the nodes along the propagation paths, suggesting a few transient propagation "scaling” regimes as a function of the nodes’ degree, but regardless of motifs such as triangles. Yet, empirical networks consist of motifs enabling the proper functioning of the system. Here, we show that basic motifs along the propagation path jointly determine the previously proposed scaling regimes of distance-limited propagation and degree-limited propagation, or even cease their existence. Our results suggest a radical departure from these scaling regimes and provide a deeper understanding of the interplay of self-dynamics, interaction dynamics, and topological properties.

Here, we study how basic undirected motifs, in particular edges and triangles, determine the response times to perturbationsas a function of the network structure and its dynamics. This allows us to identify genuine scaling regimes jointly arising from nodes' degree and motifs. Our framework of response dynamics to perturbations on networks conceptually links the small and large scale topology of a network with the spatiotemporal spreading induced by a single small perturbation. In particular, we identify and predict the impact and interplay of propagation paths, degree and motif distributions, and interaction dynamics.

Model
We characterize the network dynamics by pairwise interacting nodes, which describes the evolution of the state variables x i of node i, where throughout the manuscript i = 1…N. Different dynamics on networks are captured by the triplet {F(x i ), H 1 (x i ), H 2 (x j )} comprising of (possibly) nonlinear functions, where F(x i ) specifies the self-dynamics governing influx, leaking dynamics, degradation or reproduction [29][30][31][32][33][34] . The terms H 1 (x i ) and H 2 (x j ) determine the the adjacent interactions of node i with its neighbors, such as infection, mutualism and competition [29][30][31][32][33][34] . The connectivity matrix A accounts for the connections.
The unperturbed system is considered in a stationary collective state, which we characterize by the set of N stable equilibria, We examine the signal propagation by studying the transient dynamics induced by a permanent perturbation Δx m on the steady state x * m of the source node m. The perturbation forces nodes to transition to the shifted states x i ð1Þ = x * i + Δx i ð1Þ. For each node i, we characterize this transition period in terms of the response time τ im , defined as the time the response ratio takes the fixed value δ i (τ im ) = η. The evolution of the collective dynamics of all states may exhibit a variety of intricate spatiotemporal patterns, which we quantify by the relationship between the response time τ im and the node degree d i , the number of node i's edges. In contrast to previous work, our framework accounts for multipath connections between source and target, which are salient features of empirical networks. Figure 1 illustrates the main mechanisms underlying the impact of motifs on response times in a protein-protein network 35 , a small-world network, and an Erdös-Rényi network 11 . Basic motifs, defined as convex regular n-gons such as edges (n = 2), triangles (n = 3), squares (n = 4), and pentagons (n = 5) are ubiquitous units of random networks, as shown in Fig. 1b, c. However, as detailed in Fig. 1c, single edges disjoint from motifs are comparably rare, leading triangles to dominate networks.
To demonstrate the dynamical role of triangles we study population dynamics on networks and quantify the relative response to a perturbation Δx i ðtÞ Δx i ð1Þ with respect to different number of triangles. As shown in Fig. 1d, target nodes i as part of edges respond faster than target nodes as part of triangles, hence impacting the response time qualitatively. As shown in Fig. 1e, not only the response time may depend on the number of triangles along a path, but also the response time may vary, depending on the local network structure, even for the same number of triangles along the path. This demonstrates the impact basic motifs may have on local response dynamics on networks. a From left to right: Schematic plot of protein-protein network (shown for N = 100 nodes from total N = 2035) 35 , small-world network, and Erdös-Rényi (ER) network with network size N = 100 and average degree 10. b Share of motifs (edges, triangles, squares, pentagons) for a network. Triangles play important roles, no matter for the networks with dense or random connections, or highly clustered networks. c Share of edges as part of basic motifs (triangles, squares, pentagons). The share is calculated by the number of edges as part of corresponding motif divided by the network size. Most edges form triangles, and for the small-world network with high clustering more and more edges form larger number of triangles, demonstrating the dominance of triangles. d Target nodes i as part of the (independent) edges (brown) respond faster than other target nodes i as part of triangles (blue), as shown by Δx i ðtÞ Δx i ð1Þ . e Histogram of average propagation time from randomly chosen sources to their randomly chosen adjacent target nodes, as a function of the number of triangles (1)(2)(3)(4)(5). Error bars indicate standard deviation. d and e are based on single ER network realizations with the linking probability p = 0.10 and the network size N = 100 (rightmost network in (a)). f Perturbation response ΔxðtÞ Δxð1Þ of nodes i, j and h of the four-node network model as shown inside the panel, with target node i and its neighbors j and h, to a perturbation on the source m. Responses of j and h differ, mainly due to their different positions relative to m, but all nodes respond non-instantaneously on the same time scale, which motivated us to derive the framework, see main text. Node j responds the slowest because it is two edges apart from source m. For (d), (e) and (f), the system is governed by population dynamics _

Local propagation
To systematically quantify the impact of basic motifs on response dynamics on networks, we formulate a general theoretical framework based on Eq. (1). This requires the determination of the responses Δx i (t) to the perturbation Δx m . We are first concerned with the quantification of the local propagation from node m to its neighbor i. For small perturbations, employment of linear response theory allows us to formulate the response dynamics, where H 0 2 ðxÞ represents the derivative dH 2 (x)/dx, which is evaluated at the initial states x * j and x * m , while J i represents the self-dynamics of the form The right-hand side of Eq. (3) holds three contributions to the response: self-dynamics, a sum specifying the adjacent interaction dynamics j ≠ m, and the response of node i directly induced by source m. The perturbed system converges for t → ∞ to the system's new collective stationary state, with responses Δx i (∞). We quantify the responses in finite time by the ratio δ i (t), Eq. (2), whose solutions are obtained by reformulating the linear response equation (3), where E im ðtÞ represents the contribution from the neighbors' (adjacent) dynamics and is defined as and the response times are determined by the solutions δ i (t = τ im ) = η.
If we assumed that the states of adjacent nodes jump instantly to their new stationary state, Δx j (τ im ) ≈ Δx j (∞), the term E im ðtÞ would vanish, E im ðtÞ ≈ 0. Previous work 15,36 that has hypothesized three distinct spatiotemporal scaling regimes has implicitly assumed this premise. However, it is important to emphasize that the premise is not always valid, especially for networks with prevalent motifs, and not even for the four-node network in Fig. 1f. To recognize this, we focus on this four-node network and its responses to a perturbation on node m, Δx i (t), Δx h (t) and Δx j (t), where nodes m, i and h form a triangle, and node j is adjacent to node i. Note that edge (i-j) is referred to as an independent edge as it is not directly connected to the source m. The convergence behaviors of the responses Δx h (t) and Δx j (t) are comparable with Δx i (t), which is conflicting to a small E im ðtÞ. We also observe that Δx h (t) and Δx j (t) exhibit qualitatively different asymptotic behaviors. Also note that the response of node j is slower than those of nodes h or i because j is two edges apart from source m, which has a stronger effect than the overall faster responses of nodes as part of independent edges compared to triangles.
Apart from triangles, other motifs may also occupy a large proportion of the networks, as shown in Fig. 1b, c. To further analyze the impact of triangles and other basic motifs, we compare the response times in small synthetic networks with and without localized signal flow disruptions (see Supplementary Material, Tables S1 & S2). We study n-gons, from triangles (n = 3) to pentagons (n = 5) finding that differences of the response time are larger for the smaller n. This means that the most significant impact on the response time can be attributed to independent edges and triangles, which prompt us to decompose networks into independent edges, that is, 2-gons such as the (i-j)-edge in the example in Fig. 1f, and triangles (3-gons).
We quantify the response times from three basic perspectives: self-dynamics, independent edges and triangles. The unperturbed system is considered in a steady state, characterized by N stable equilibria x * i . The relationship between the intrinsic dynamics of node i and its adjacent dynamics in the steady state follows immediately from Eq. (1) and _ x i = 0, Averaging allows us to compute the contribution of adjacent dynamics in the mean-field which is then used to simplify the relationship (7) as in which the degree serves as a coupling constant of the intrinsic dynamics of node i and its adjacent connections. In the unperturbed steady system, the equilibria x * i can be expressed as the inverse function of Rðx * i Þ with respect to the degree d i , Combining Eqs. (5) and (9), we derive the response time τ im of node i as a function of its degree d i , as where , the quantity Q i is a function of degree d i , while the mean-field quantity Q im is independent of the degree. In the Supplementary Material, the Hahn expansion of Q i leads to its leading power, Q i~d with the scaling exponent θ Q = Π Q (∞) + 1. The exponent θ Q is determined by the intrinsic dynamics but is independent of d i . In the large-degree limit, for θ Q < 0 the contributions to τ im from adjacent nodes vanish such that the response time becomes independent of m and is well approximated by This result coincides with existing literature and mechanistically explains the validity of previous theoretical results for a number of dynamics in the large-degree limit 15,36although the adjacent interactions as quantified by E im are not considered.
The self-dynamics term J i can be expanded as a function of the degree d i , as J i~d θ J i , where the scaling exponent θ J is the leading power of the Hahn expansion. In doing so, we find that the response time τ i exhibits the scaling relation where the scaling exponent θ J is determined by the intrinsic dynamics but independent of adjacent connections. The scaling relationship (13) highlights the contribution of both the structural features and system dynamics on the response time, yet it is a disentanglement of self-dynamics J i and degree d i . In that way, this finding is in agreement with previous work on three distinctive dynamic regimes 15 . In contrast, for θ Q > 0, the contributions from adjacent nodes may be substantial, even for large degree. In this case, both the self-dynamics and the adjacent dynamics contribute to the response time, remarkably exhibiting a scaling relation where the scaling is affected by the adjacent dynamics, through the exponent θ Q . It is important to relate the exponent θ Q to prototypical network dynamics. As derived in the Supplementary Material, we find that regulatory, human, mutualistic, biochemical and epidemics dynamics are characterized by θ Q < 0, whereas population and inhibitory neuronal dynamics show θ Q > 0, as shown in Table 1. As the model parameters are assumed to be non-negative but otherwise arbitrary, Table 1 characterizes a substantial range of dynamical systems. The theoretical derivation, Eq. (10), exactly predicts the response time τ i . The scaling relationships, Eqs. (13) and (14), characterize the interplay between network topology, the self-dynamics and its adjacent dynamics in the asymptotic regime d i → ∞. We find that both the theoretical derivation and the scaling predictions are in good agreement with simulations for regulatory and population dynamics, as supported by Fig. 2a-e. For regulatory dynamics, characterized by θ Q < 0, the response time τ i is Table 1 | Scaling exponents θ J and θ Q for prototypical dynamical models  Fig. 2b. For population dynamics, characterized by θ Q > 0, we find three distinctive dynamic regimes. In the degree-limited regime with θ J > 0 in Fig. 2c and the distance-limited regime with θ J = 0 in Fig. 2d, the term E im is negligible due to the degree limitation, and scaling is only determined by the self-dynamics Eq. (13). However, in the composite dynamic regime (θ J < 0), τ i is determined by both the self-dynamics and the adjacent dynamics as predicted by Eq. (14) and supported by Fig. 2e. While existing literature disregarded the effects from adjacent dynamics and predicts onlyθ = θ J 15 , leading to inaccurate scaling as shown in Fig. 2e, our prediction, θ = θ J − θ Q , is asymptotically exact. Taken together, we observe scaling of the response time that depend on both the self-dynamics and the adjacent dynamics, even for local tree-like networks. But how relevant is this composite dynamic regime?
For a triangle-dominated topology as shown in Fig. 2f, we derive the response time τ i from Eqs. (5) and (9) as where C im can be approximated by the product of the degreeindependent mean-field term Q im and d i , C im ≈d i Q im . Note that f is a constant that depends on the system dynamics, but not on d i . Now, for f ≪ 1 and a large degree, C im drops out in Eq. (15). In this case, triangles have a negligible impact on the scaling and the scaling relation is well approximated by τ i~d θ J i , as shown in Fig. 2g. For f ≈ 1, the presence of triangles, as shown in Fig. 2h, results in the extra factor d i resulting from C im such that As shown in Fig. 2f-h for regulatory and population dynamics, both predictions, Eqs. (15) and (16) are well confirmed by our computer simulations. This establishes that basic motifs may crucially determine scaling, even in simple networks.
In the Supplementary Material, we derive the respective explicit solutions of the response time, not only for the asymptotic scaling regime but also for the regime of small degrees, which allows us to fully characterize the scaling regimes. In addition to regulatory and population dynamics, we develop the system dynamics framework for human 29 , epidemic 30 , mutualistic 31 , biochemical and inhibitory dynamics [32][33][34]37,38 (see Table 1). Our results consistently quantify the impact of basic motifs on local response dynamics on networks.

Global propagation
Thus far, we have established a versatile theoretical framework for local signal propagation. To formulate a framework for global signal propagation, it is helpful to examine a source-centric representation, which allows to theoretically track all possible propagation paths from the perturbed central node to all distant nodes. Specifically, we impose a layer-to-layer topology with the source node in the center and study the interplay of the intrinsic and adjacent dynamics, employing our developed framework for the local propagation as a building block. Denote by T(m → i k ) the propagation time from node m to the i k -th node through a pathway, specifying when the signal response ratio attains the threshold η. The signal propagation is expressed layer-wise in terms of T(m → i k ), which we compute recursively, where the subscript k of index i k stands for the distance (the number of edges along the path) to the perturbed node m. Based on the Gauss Iterative Method combined with the generalization of Taylor expansion 39 , we solve Eq. (17), yielding with vector g T = gð1Þ gð2Þ . . . gðkÞ ð Þ , whose components are monotonic decreasing functions that approach 1 as k → ∞. The components g(h) of g are approximately proportional to a product, accounting for the accumulated dynamical effects from independent edges of layers higher than h. Note that the components may vary across different system dynamics, but are of the same order for a given dynamics (see Supplementary Material). The vector D depends on the degree sequence d i k and the parameters characterizing the adjacent dynamics, where C 1 and C 2 are constants, and the exponent θ J characterizes the self-dynamics, while θ Q depends also on the adjacent dynamics, see Table 1. The scalar product (18) is the summation of the propagation times through the layers, where the vector g weights the degree sequence along the pathway as held in D. Because the components of g are given as product (19), degree fluctuations do not average out as it would for an additive structure. This is why and how the degree sequence affects the global propagation time, T(m → i k ). Taken together, equation (18) predicts that not only the average degree, as assumed in previous literature 15 , but also the variance and even the degree sequence along the propagation path may crucially determine the response times to a perturbation. To further illustrate this, we study how signal propagation is impacted by the average degree and degree variance of the nodes forming a chain. Figure 3a, c show for regulatory dynamics with θ J > 0 that the propagation time increases with increasing mean degree, and decreases with increasing standard variation of the randomly chosen degree sequence. In contrast, for θ J < 0, the propagation time decreases with the mean degree and the standard variation, as shown in Fig. 3b, d.
As we expect the global signal propagation T(m → i k ) to be sensitive to the degree sequence along the propagation path, in the Supplementary Material, we systematically study degree sequences with large standard variation, which leads us to two crucial observations. First, for θ J > 0, the propagation time is dominated by the largest node degree maxðd θ J i Þ along the chain, which is consistent with the numerics as shown in Fig. 3e. Second, for regulatory dynamics with θ J < 0, the propagation time is dominated by the smallest degree minðd θ J i Þ, which is supported by Fig. 3f. To study the impact basic motifs have on empirical networks, we analyzed a protein-protein network 35 , where we control the clustering coefficient by edge rewiring. The central node is the source node inducing a perturbation. Nodes in the same layer (with the same radial distance) have the same shortest path length to the source, as shown in Fig. 4a. Brown nodes indicate the arrival of the perturbation in (an arbitrarily) fixed time, for fixed η. For this dynamics, triangles inhibit propagation. The average propagation time increases as a function of number of triangles and layers, as shown in Fig. 4b. Similarly, we find that signal propagation is increasingly inhibited with increasing clustering coefficient, as shown in Fig. 4c. For an extensive analysis of this empirical network, together with a detailed investigation for a number of prototypical networked systems, whose arbitrary parameter ranges are expected to characterize a wide range of empirical networked dynamical systems, we refer to the Supplementary Material.

Discussion
Triangles and other loops have always limited the understanding of the interplay of function and structure in networks. We have developed analytical tools that allow to capture the impact of simple undirected motifs on the system dynamics. Our developed framework not only helps disentangle joint effects but provides a deeper understanding of the interplay of self-dynamics, interaction dynamics, and topological properties. Our analysis suggests a radical departure from the previously proposed concepts of distance-limited propagation and degree-limited propagation. In distance-limited propagation the response time scaling is said to be dominated by the propagation path length but not by the edge density along the path. Vice versa, for degree-limited propagation the response time scaling is said to be dominated by the mean degree, but not the propagation path length. We have demonstrated here by independent methods that when the propagation is drastically slowed, or accelerated, that may not necessarily result from edges or hubs but from cycles, in particular triangles. Our analysis is based on a network decomposition into independent edges and edges as part of motifs. The developed framework predicts genuine scaling exponents, no matter if the propagation dynamics is dominated by hubs, by the path length, or by basic motifs. For paths with large average degree, as abundant in social and other empirical networks, the prediction of propagation time using asymptotic scaling as proposed by existing literature may be orders of magnitude off, as the scaling exponent may fall in an unrelated universality class. We have overcome this inconsistency by introducing two topologyindependent exponents that quantify the universality class of the local response dynamics on networks.
Network motifs are abundant in synthetic and empirical networks and systematically impact the response dynamics to perturbations. We have provided a versatile toolbox that may not only help understanding response dynamics on networks but also provide the mathematical building blocks for extensions such as genuine universality classes for dynamics on directed multiplex networks. Yet, it is a Signal propagation snapshots of protein-protein networks 35 . Central node is the source node inducing a perturbation. Nodes in the same layer (same radial distance) have same shortest path length to the source. Brown nodes indicate the arrival of the perturbation in (an arbitrarily) fixed time. b Average propagation time as a function of number of triangles and layers in (a). For a given layer, triangles slow propagation. Average propagation time increases with both number of triangles and layers. c Signal propagation in rewired networks with varying the clustering coefficient by edge rewiring. For all three networks, high clustering and triangle density slows the spreading of the perturbation across layers. Panels for regulatory dynamics, _ x i ðtÞ = À Bx i a + α ∑ N j = 1 A ij important to note that the developed analytical tools are built on linear response theory to quantify the dynamics in the steady state as a linear response to a small permanent perturbation of a single unit. Large dynamic perturbations, for instance, may force the system to transition into steady states not predicted by linear response theory. Theoretical extensions covering this subject remain challenging but deserve future attention.

Data availability
The data in this study is freely accessible at https://github.com/ QitongHu2000/Impact-of-motifs-data.